#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkSliceBySliceImageFilter.h"
#include "itkBinaryThinningImageFilter.h"
#include "itkCommand.h"
#include "itkSimpleFilterWatcher.h"
#include "itkLabelObject.h"
#include "itkLabelMap.h"
#include "itkBinaryFillholeImageFilter.h"

/* convexHull.hdr and walloutput.hdr are both region images;

./PhaseSpliter 20 3 0 0.9722222089767456 7.0 ./convexHull.hdr;
./PhaseSpliter 20 3 0 0.9722222089767456 7.0 ./walloutput.hdr;

./MidMyocardiumConstruction convexHull_8.hdr walloutput_8.hdr midMyocardium_8Contour.hdr;
./MidMyocardiumConstruction convexHull_20.hdr walloutput_20.hdr midMyocardium_20Contour.hdr;

./fillVolume midMyocardium_8Contour.hdr midMyocardium_8.hdr 0 1;
./fillVolume midMyocardium_20Contour.hdr midMyocardium_20.hdr 0 1;


./SliceInterpolator convexHull_8.hdr 0.9722222089767456 7.0 convexHull_8Interpolated.hdr;
./SliceInterpolator midMyocardium_8.hdr 0.9722222089767456 7.0 midMyocardium_8Interpolated.hdr;
./SliceInterpolator walloutput_8.hdr 0.9722222089767456 7.0 walloutput_8Interpolated.hdr;

./SliceInterpolator convexHull_20.hdr 0.9722222089767456 7.0 convexHull_20Interpolated.hdr;
./SliceInterpolator midMyocardium_20.hdr 0.9722222089767456 7.0 midMyocardium_20Interpolated.hdr;
./SliceInterpolator walloutput_20.hdr 0.9722222089767456 7.0 walloutput_20Interpolated.hdr;

./GenParaMesh  --label 1  --iter 200  --inVolName ./convexHull_8Interpolated.hdr   --outParaName ./convexHull_para8.vtk   --outSurfName ./convexHull_surf8.vtk;
./GenParaMesh  --label 1  --iter 200  --inVolName ./midMyocardium_8Interpolated.hdr   --outParaName ./midMyocardium_para8.vtk   --outSurfName ./midMyocardium_surf8.vtk;
./GenParaMesh  --label 1  --iter 200  --inVolName ./walloutput_8Interpolated.hdr   --outParaName ./walloutput_para8.vtk   --outSurfName ./walloutput_surf8.vtk;

./GenParaMesh  --label 1  --iter 200  --inVolName ./convexHull_20Interpolated.hdr   --outParaName ./convexHull_para20.vtk   --outSurfName ./convexHull_surf20.vtk;
./GenParaMesh  --label 1  --iter 200  --inVolName ./midMyocardium_20Interpolated.hdr   --outParaName ./midMyocardium_para20.vtk   --outSurfName ./midMyocardium_surf20.vtk;
./GenParaMesh  --label 1  --iter 200  --inVolName ./walloutput_20Interpolated.hdr   --outParaName ./walloutput_para20.vtk   --outSurfName ./walloutput_surf20.vtk;

./ParaToSPHARM  --degree 10  --inParaFile ./convexHull_para8.vtk   --inSurfFile ./convexHull_surf8.vtk  --outDir ./ --outSurfFile convexHull_SPHARM8.vtk;
./ParaToSPHARM  --degree 10  --inParaFile ./midMyocardium_para8.vtk   --inSurfFile ./midMyocardium_surf8.vtk  --outDir ./ --outSurfFile midMyocardium_SPHARM8.vtk;
./ParaToSPHARM  --degree 10  --inParaFile ./walloutput_para8.vtk   --inSurfFile ./walloutput_surf8.vtk  --outDir ./ --outSurfFile walloutput_SPHARM8.vtk;

./ParaToSPHARM  --degree 10  --inParaFile ./convexHull_para20.vtk   --inSurfFile ./convexHull_surf20.vtk  --outDir ./ --outSurfFile convexHull_SPHARM20.vtk;
./ParaToSPHARM  --degree 10  --inParaFile ./midMyocardium_para20.vtk   --inSurfFile ./midMyocardium_surf20.vtk  --outDir ./ --outSurfFile midMyocardium_SPHARM20.vtk;
./ParaToSPHARM  --degree 10  --inParaFile ./walloutput_para20.vtk   --inSurfFile ./walloutput_surf20.vtk  --outDir ./ --outSurfFile walloutput_SPHARM20.vtk;
./MidMyocardiumConstruction convexHull_20.img walloutput_20.img midMyocardium_20Contour.img
*/

int main(int argc, char * argv[])
{

  if( argc != 4 )
    {
    std::cerr << "usage: " << argv[0] << " endocardium epicardium  midMyocardium" << std::endl;
    std::cerr << "argc =" << argc << std::endl;
    exit(1);
    }

  const int dim = 3;
  
  typedef itk::Image< unsigned char, dim > ImageType;

  typedef itk::ImageFileReader< ImageType > ReaderType;
  ReaderType::Pointer reader1 = ReaderType::New();
  reader1->SetFileName( argv[1] );
  ReaderType::Pointer reader2 = ReaderType::New();
  reader2->SetFileName( argv[2] );

 
  typedef itk::BinaryThresholdImageFilter <ImageType, ImageType> BinaryThresholdImageFilterType;
 
  BinaryThresholdImageFilterType::Pointer thresholdFilter1 = BinaryThresholdImageFilterType::New();
  thresholdFilter1->SetInput(reader1->GetOutput());
  thresholdFilter1->SetLowerThreshold(1);
  thresholdFilter1->SetUpperThreshold(100);
  thresholdFilter1->SetInsideValue(1);
  thresholdFilter1->SetOutsideValue(0);

  BinaryThresholdImageFilterType::Pointer thresholdFilter2 = BinaryThresholdImageFilterType::New();
  thresholdFilter2->SetInput(reader2->GetOutput());
  thresholdFilter2->SetLowerThreshold(1);
  thresholdFilter2->SetUpperThreshold(100);
  thresholdFilter2->SetInsideValue(1);
  thresholdFilter2->SetOutsideValue(0);


  typedef itk::SubtractImageFilter <ImageType, ImageType > SubtractImageFilterType; 
  SubtractImageFilterType::Pointer subtractFilter = SubtractImageFilterType::New();
  subtractFilter->SetInput1(thresholdFilter1->GetOutput());
  subtractFilter->SetInput2(thresholdFilter2->GetOutput());
  subtractFilter->Update();


  typedef itk::SliceBySliceImageFilter< ImageType, ImageType > SliceBySliceImageFilterType;
  SliceBySliceImageFilterType::Pointer sliceBySliceFilter = SliceBySliceImageFilterType::New();
  sliceBySliceFilter->SetInput( subtractFilter->GetOutput() );


  typedef itk::BinaryThinningImageFilter < SliceBySliceImageFilterType::InternalInputImageType, SliceBySliceImageFilterType::InternalInputImageType >  BinaryThinningImageFilterType;
  BinaryThinningImageFilterType::Pointer  binaryThinningFilter = BinaryThinningImageFilterType::New();

  sliceBySliceFilter->SetFilter( binaryThinningFilter );

  typedef itk::ImageFileWriter< ImageType > WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput( sliceBySliceFilter->GetOutput() );
  writer->SetFileName( argv[3] );
  writer->Update();
  return 0;
}


